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Abstract 


Element free Galerkin (EFG) method is a computational tool in solid mechanics 
problem for solving partial differential equations based on a scattered set of nodal 
points. This method belongs to the class of meshless methods, which is particularly 
useful for problems involving crack propagation due to the absence of any 
predefined element connectivity, and the burdensome remeshing required as in 
finite element method (FEM). This method is based on the use of moving least 
square Interpolants with a Galerkin method. In the present problem, the 
implementation of EFG method for problem of fracture and quasi static crack 
growth is described. Numerical results show that accurate stress intensity factor 
can be obtained without any enrichment of displacement field near the crack tip. 
The stress singularity and the crack growth can be easily modeled. In this crack 
propagation problem in plane elasticity is considered. Lastly, the problem of 
propagation of multiple, arbitrarily placed, crack is considered. 
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CHAPTER 1 


INTRODUCTION 


1.1 INTRODUCTION 

For the many structures, crack propagation is an important failure mechanism 
requiring accurate numerical models to implement simulations essential for failure 
prediction. To perform numerical simulation, the computational methods must be 
applied to determine fracture response and reliability of cracked structures. A 
current popular method is finite element (FEM), which has been extensively used 
for fracture analysis of cracks. Although a significant amount of the research in 
FEM is useful, the method has serious limitations in solving solid mechanics 
problems characterized by a continuous change in geometry of the domain under 
analysis. Crack propagation is a prime example in which the use of FEM requires a 
large number of remeshings of the finite element model to represent arbitrary and 
complex paths. The underlying structures of FEM and similar methods, which rely 
on a mesh, is quite cumbersome in treating cracks, that is not coincident with the 
original mesh geometry. Consequently, the only viable option for dealing with 
moving cracks model evolution so that mesh lines remain coincident with cracks 
throughout the analysis. This creates numerical difficulties, often leading to 
degradation of solution accuracy, complexity in computer programming, and a 
computationally intensive environment. 

The principal attraction of meshless methods is possibility of simplifying 
adaptively and problems with moving boundaries and discontinuities, such as 
phase changes or cracks. In crack growth problems, for example, nodes can be 
added around a crack tip to capture the stress intensity factor with the desired 
accuracy; this nodal refinement can be moved with a propagating crack through a 
background arrangement of nodes associated with the global geometry. Adaptive 
meshing for a large variety of problems, including linear and non-linear stress 
analysis, can be effectively treated by these methods in a simple manner. 
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1.2 LITERATURE SURVEY 


In the field of computational mechanics, extensive research has been done in the 
effective use of finite element method for problems in solid mechanics (fracture 
mechanics). Effort has been primarily focused on the modification of the 
approximation, incorporating the singularity present in the domain. The element 
free Galerkin method is a recently developed computational tool, which has been 
extensively for computing arbitrary crack paths. This class of methods is often 
called meshless, gridless or particle methods because of the absence of any 
predefined nodal connectivity. 

In the past decades, for the solution of growing crack problems finite 
element methods with remeshing and boundary element methods is used but each 
has substantial disadvantages, which preclude their application to many important 
engineering problems. In the first one that involves a considerable amount of 
engineering effort and can result in substantial errors, it is almost impossible to 
automatically remesh finite elements about arbitrarily growing crack. Boundary 
element methods avoid a large part of the remeshing associated with finite element 
methods because only the surface of the crack, or in two dimensions, the line 
representation of the crack, needs to be remeshed. However, since boundary 
element methods require a Green function for the underlying partial differential 
equations, they are quite limited in the scope of problems they can attack. 

Meshless method has been in existence since the late 1970’s. Lucy (1977) 
introduced a particle method called smoothed particle hydrodynamics (SPH) for 
modeling astrophysical phenomena and in the 1980’s Monaghan (1982) used this 
method for problems without boundaries, such as rotating stars and dust clouds 
and to solve the solid mechanics problems and further extended to solve the solid 
mechanics problems. In recent years, a class of meshfree or meshless methods, 
such as smooth particle hydrodynamics, diffuse element method, element free 
Galerkin method, hp clouds, partition of unity, reproducing kernel particle method, 
meshless local Petrov-Galerkin method and local boundary integral equation 


2 



method appear to demonstrate significant potential for the moving boundary 
problem like growing cracks. Fundamental to all mesh free methods, a structured 
mesh is not used, since only a scattered set of nodal points is required in the 
domain of interest. 

Meshless methods based on partitions on unity seem to provide an efficient 
vehicle for performing hp adaptivity. These methods include hp -clouds and the 
partition of unity finite element method (PUFEM). The hp -clouds methods begins 
with a partition of unity based on moving least squares and enhances the 
polynomial order of the approximation through an extrinsic basis which can be 
added locally to nodes. 

Other meshless methods, which have evolved, are the particle in cell (PIC) 
method, the generalized finite difference method and the finite point method. 
Belytschko (1996) provide a comprehensive review of the state-of-the-art in 
meshless methods and details the relationship between several of the key 
methods. 

The essential idea of EFG methods is the use of moving least squares 
interpolants. Such interpolants were first described in a limited form by Shepard 
(1968): a more complete description can be found in the article by Lancaster and 
Salkauskas (1981). These interpolants were evidently independently discovered 
later who were the first to sue these in Galerkin methods: they called them diffuse 
elements and used these methods or the solution of heat conduction problems. 

A formulation has been presented by Belytschko (1995) for consistently 
coupling EFG and FE by blending the approximations in a transition element. This 
allows the speed and simplicity of finite elements to be exploited where possible 
while allowing EFG to be used in regions where a meshless method is appropriate. 
But our problem of crack propagation for arbitrary paths is inherently difficult and 
also compromises the salient feature of EFGM. Jasiuk (1995) modeled fracture and 
crack growth with finite elements by deleting which met a yield criterion. This 
approach, which is not based on fracture mechanics, requires a very fine mesh to 
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get an acceptable representation of a crack. Other techniques for modeling crack 
growth include spring network models in which the material is represented by a 
network of spring and crack propagation is simulated by breaking spring Schlangen 
and Mier (1992). 

The most obvious way to handle crack propagation is by remeshing the 
geometry. Swenson and Ingraffea (1988) presented a local remeshing technique in 
which elements ahead of the crack tip in the propagation direction were removed 
and the crack was extended. This area was triangulated to create a new local crack 
tip mesh. This method has the advantage that a vast amount of existing finite 
element technology can be used for support. There are some drawbacks to this 
method such as difficulties if the crack step size is too small, then remeshing is 
unwarranted. Complex geometries and interacting crack tip are also potential 
difficulties as is rerunning a step with a smaller crack increment. 

1.3 Thesis organization 

This thesis is focus on the element free Galerkin (EFG) methods for computational 
fracture mechanics and is organized as follows: 

In chapter 2, we briefly review of linear elastic fracture mechanics and introduce 
the J integral for numerically computing stress intensity factors. Quasi-static 
fracture crack is discussed and algorithm for their simulation by element free 
Galerkin (EFG) are introduced. 

In chapter 3, the moving least square methodology is reviewed and used to 
construct discrete EFG approximation. The elastic static boundary value problem is 
presented along with its associated weak form. Approximation of the solution with 
EFG is presented and topics of nodal domains of influence and integration of weak 
form are discussed. It should be noted that the definition of a meshless method is 
one in which no predefined element connectivity exists for determining the 
approximant. Some confusion invarbly arises when the meshless approximant is 
used in Galerkin method. Integration of weak form is performed by Gauss 
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quardrature, which requires some sort of integration cells. Although this detracts 
from the meshless nature of the method the background cell structure by no means 
destroys it. Smoothing of EFG approximation near nonconvex boundaries is 
presented. Without smoothing, EFG approximation near nonconvex boundaries 
such as crack tips will contain discontinuities, which extend into the domain. These 
discontinuities arise due to a sampling point grazes the boundary. The diffraction 
method smoothes EFG approximation by wrapping the nodal support a short 
distance around the point at which the discontinuity would begin. 

In chapter 4, several example problems are presented to illustrate the 
effectiveness of EFG for solving problem of arbitrary crack growth. 

Finally in chapter 5, the work is summarized and the relevant conclusions and the 
scope of the future work were presented. 
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CHAPTER 2 


LINEAR ELASTIC FRACURE MECHANICS 

2.1 Parameterization of fracture field 

The presence of a crack has a considerable effect on the stress field and load 
bearing capacity of a body. For cracks in linear elastic media, the stress fields are 
well known and can be expressed in terms of parameters, which include the effects 
of crack length, applied loads, and other conditions such as boundary effects. The 
two parameters, which will be discussed, are the stress intensity factor and the 
energy release rate. The stress intensity factor is a local crack tip parameter, which 
describes the stress field at a crack tip, while the energy release rate is a global 
parameter, which gives the change in elastic energy for an infinitesimal crack 
advance. 

Fracture modes 

The displacement field around a crack can be categorized by the three 
modes shown in Fig 2.1. Mode I is called the opening mode and is charactized by 
displacements of the crack normal to the plane of the crack. Mode II is called the 
shearing or sliding mode and is characterized by in-plane displacements of the 
crack faces parallel to the crack surface. Mode III is called the tearing or antiplane 
mode and is characterized by shear deformation out of the plane of the crack. 

In this report, only two dimensional plane problems are presented and 
hence mode I and mode II fracture are considered. 

Stress intensity factors 

The stress at the tip of a linear elastic crack is 1 l4r singular, where r is the 
distance from the crack tip. The stress field near the crack tip is dominated by the 
first term in the asymptotic. For a mode I crack (displacement perpendicular to the 
crack faces), the stresses can be written as 
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Opening Mode Shearing Mode 



Tearing Mode 


Fig. 2.1 Schematic drawing of the three fracture modes 
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The variables r and 0 are the radius from the crack tip and the angle 
measured from the line ahead of the crack, respectively (see Fig. 2.2) is the 
Kolosov constant k written in terms of Poisson’s ratio v . 


k = 


3-4 v plane strain 
3 -u 

plane stress 

1 + o 




Fig. 2.2 Local coordinate system at crack tip. 


The mode I and mode II stress intensity factors, K, and K„ are 

constants, which depend on the crack length, specimen geometry and applied 
loading. It is significant to note that determination of these constants completely 
determines the asymptotic stress and displacement fields around the crack tip. 

In case where the loading and geometry give a combination of opening 
(mode I) and shearing (mode II) displacements, the stress field will be a 
combination of previous equations. In polar coordinates, the crack tip stresses are 
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2.2 Computation of fracture parameters 

A great deal of research has gone into methods of determining stress intensity 
factor for various geometries and loadings. The developments of numerical 
methods such as finite elements and boundary elements have facilitated the 
calculation of crack tip stress fields and many methods have been developed for 
inferring stress intensity factor (e.g., displacement and stress extrapolation, 
displacement and stress correlation, enriched elements, J integral). The J integral 
is popular because of its robustness as well as adaptability to mixed mode 
problems. 

J integral 

The J integral is given by 

. . . .(2,14) 

r L ® x \ J 


where W = j ads is the strain energy density, a is the Cauchy stress, and n is 

0 

the outward unit normal to an arbitrary contour enclosing the crack tip. The path 
independence property allows the J integral to be evaluated using far field 
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information, which is generally more accurate, that near-tip results for numerical 
methods. Path independence is easily shown by considering the closed contour 

T = T, +C, +F, +C_ Around the tip of a straight crack with m as the normal to the 
closed contour. Along the entire contour previous can be written 
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where the prime is used to indicate that J is written for a closed contour. Along the 
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last integral equation is zero. Reversing the direction of the normal m on the inside 
contour T, to be the outward normal n r leads to J r =J r . 
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Domain form of J integral 

The contour integral equation can be converted into equivalent form to 
enhance the usefulness of the J integral in certain cases. The evaluation of contour 
integrals is difficult with finite elements because the stress field is discontinuous 
and performing an area integral is more convenient. 


Consider a closed contour with weight function q included in the integeral and 
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Since q-0 on r 2 and the integer and vanishes on C+ and C., the closed 
contour integral is reduces to a closed integral on r, 
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Fig. 2.3 Path Independent closed contour about the tip of a crack 

Applying the divergence theorem to obtain the area integral 
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For linear elasticity, the first term is vanishes. Comparing both equation and 
reversing the normal m in the equation to be outward normal n r leads to 
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The domain form of the J integral is then 
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Interaction integral 

The J integral defined in equation is equal to the energy release rate for 
elastic materials. For a general mixed-mode case 

J = G=±( K f +K ;) . . . .(2.21) 


where 



plane stress 
plane strain 


The individual stress intensity factors K, and K„ are needed to compute 
the direction and growth of fracture crack. 


Considering two independent equilibrium states of a cracked body, where 
state 1 is defined here to be the actual state for the given boundary conditions 
while state 2 is defined to be an auxiliary state, which will be given later. The J 
integral for the two superposed states is 
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Obtaining the individual stress intensity factors from the actual field requires 
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that the auxiliary field be chosen judiciously. For the case, k/ 2) = = 0 earlier 

equations reduce to = — K. (l) 

E 1 


Similarly, the mode II stress intensity factor can be found with 

k/ 2 ^ = 0,K n ^ = 1 equation is reduces to M ^ = — K,/ 2) 

E 


Thus, the contour integral can be converted to an equivalent domain integral 
by Gauss’s theorem and leads to the equation 
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It should be noted that interaction integral as described is only valid for 
elastic materials because of the auxiliary field chosen. However, the stress 
intensity factor as a fracture parameter is only valid when crack tip plasticity is 
negligible and so this limitation on the interaction integral should not be of concern. 

2.3 Crack growth direction 

Modeling crack propagation requires a law for the direction of crack growth. The 
most commonly used laws are 

1 . The maximum principal stress criterion 

2. The maximum energy release rate criterion 

3. The minimum strain energy density criterion 

The first two criteria have been identical for the electrostatic fracture, 
although they are not necessarily the same for dynamic fracture. The minimum 
strain energy density criterion assumes that the strain energy density near crack tip 
is a quadratic function of mode I and mode II stress intensity factors. Crack 
extension is assumed to occur in the direction, which minimizes the strain energy 
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density. In this report, the maximum principal stress criterion is used for computing 
the direction of crack' extension. 

The maximum principal stress criterion states that a crack will grow in a 
direction perpendicular to the maximum hoop stress, a gg . For a general mixed- 
mode state of stress, the stresses at the crack tip are given by the earlier equation. 
The direction of maximum hoop stress is the orientation in which the shear stress 
vanishes 


cr 
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V2 nr 2 


^K,sin| + ^K„(3cos0-l) 


. . . (2.33) 
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The solution cos- = 0 y ields an angled = ±180°, which is ignored because it 

is irrelevant. The other solution is obtained by setting the term in the square 
brackets to zero. This leads to the condition 

K, sin0 m + K I] (3cos#„ l -l)=0 - (2.34) 

which must be solved for the angle of crack propagation, 6 m . 


For this angle, a og is principal stress field. This means that a 09 may be 
written in terms of an equivalent mode I stress intensity factor instead of the mixed 
mode stress intensity factors, which leads to 
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The equivalent mode I stress intensity factor, K Ie(/ , provides a single 
measure of the mixed-mode stress field and is used in crack growth laws for 
computing the crack growth rate mixed-mode loading. 
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Multiple cracks 


The crack growth methodology presented in the previous section requires 
some modification if more than one crack is present and growing. The stress field 
at the tip of each crack will be different and thus some cracks will grow faster than 
others and some may not grow at all. 

Algorithm used for growing multiple quasi-static crack 

1 .compute the stress field 

2.compute stress intensity factors for each crack 

# Compute incipient crack growth angle 9 m using maximum principal stress 
criterion 

#Compute equivalent mode I stress intensity factor K |C(/ 


3.find crack with highest K lw/ , i.e. fastest growing crack 
Set A a = Aa imx for this crack 


4. compute the crack increment for each based on A a = A a max 


f 


v 


K* 


\ 


J 


2 


Note: if K la/ <K C for a crack, then A a - Ofor that crack 
5. increment each crack 

6. if another step of crack growth is desired, go to step 1 
7.end analysis 
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Quasi-static crack propagation 

Quasi-static fracture refers to crack propagation at or above the fracture 
toughness. No accepted standard for computing the crack growth rate exists in this 
condition, and the crack is assumed to propagate as long as the stress intensity 
factor is greater than the fracture toughness, i.e. whenever K > K /c . is satisfied. 

Crack growth may be classified as stable and unstable based on the change 
in stress intensity factor with growth. In stable fracture, the stress intensity factor 
decreases as the crack propagates and crack will eventually arrest when it falls 
below the fracture toughness. Examples of stable fracture are cracked specimens 
loaded by prescribed displacement and cracks in an overall compressive field. In 
unstable fracture, the stress intensity factor increases with crack length and 
propagating crack never arrests. 

Fracture stability is easily understood by considering a compact tension 
specimen with a prescribed displacement and prescribed load. For the prescribed 
load, the stress intensity factor continues to increase as the crack extends, 
indicating unstable fracture. For the prescribed displacement, the stress intensity 
factor decrease as the crack grows and eventually falls below the fracture 
toughness. The crack arrests at this point. 

2.4 Closure 

In this chapter briefly review of linear elastic fracture mechanics and introduce the 
J integral for numerically computing stress intensity factors. Quasi-static fracture 
crack is discussed and algorithm for their simulation by element free Galerkin 
(EFG) are introduced. 
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CHAPTER 3 


ELEMENT FREE GALERKIN METHOD 

3.1 Introduction 

The element-free Galerkin (EFG) method is a meshless method because only a set 
of nodes is required to generate the discrete equations. The connectivity between 
the nodes and the approximation functions are completely constructed by the 
method. 

The EFG method employs moving least-square (MLS) approximants to 
approximate the function u(x ) with u''(x) Belytschko, (1998). The approximation of 
z/'(x) the function z/(x)at any point x in the problem domain Q is written 

7 JT 

u (x) = p (x)a(x) . These approximants are constructed from three 

components: a weight function of compact support associated with each node; a 
basis, usually consisting of a polynomial; and a set of coefficients that depend on 
position. The weight function is nonzero only over a small sub domain around a 
node, which is called its support. The support of the weight function defines a 
node’s domain of influence, which is the sub domain over which a particular node 
contributes to the approximation. The overlap of nodal domains of influence defines 
the nodal connectivity. 

One attractive property of MLS approximants is that their continuity is 
related to the continuity of the weight function, if the continuity of the basis is 
greater than the continuity of the weight function, the resulting approximation will 
inherit the continuity of the weight function. Therefore, a low order polynomial 
basis, e.g., a linear basis, may be used to generate a higher order continuous 
approximations by choosing an appropriate weight function. Thus, post processing 
stress and strain results to generate smooth fields that is required for finite element 
methods is unnecessary in EFG. 
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Although EFG is considered meshless with reference to shape function 
construction or function approximation, a mesh will be required for solving partial 
differential equations (PDE s) by Galerkin approximation procedure. In order to 
compute the integrals in the weak form, either a regular background mesh or a 
background ceil structure is used. 

3.2 Meshless approximations by MLS 

The approximation of u h (x ) the function w(x)at any point x in the problem 
domain Q. is written 


u h (x) = p r (x)a(x) (3.1) 

where p(x ) is a basis (usually polynomial) and a(x) are unknown coefficients: 
Some complete polynomial basis is 


/? r (x) = [l,x] 

linear. 

(3-2 a) 

p‘ (x) = [l,x,x 2 ] 

quadratic. 

(3-2 b) 

//(x) = [l,X,y] 

linear. 

(3.3 c) 

p' (x) = h,x,y,x 2 ,xy,y 2 ] 

quadratic. 

(3.2 d) 


Sometime some other function can be added to the basis in situations where 
it is desirable to enrich the solution. 

To find field variable by earlier equation, it is necessary to determine the 
coefficients a(x ) . The moving least squares (MLS) methodology [Lancaster and 
Salkauskas, 1981] is an effective technique for approximating a function using a 
set of scattered data. Given a set of nodes with coordinates x, at which the field 
variable u , is known, weighted norm L 2 can be written 
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. (3.3) 


J(x) = J w(x - x , )[p T (x, )a{x) - u, ? 

/ = 1 

where w{x- X[ ) is the weight function of node I at point x, and n is the number of 
neighbours of point x, i.e., nodes with w(x-x / )>0 and u, refers to the nodal 
parameter of u at x = x, Note that u h (x,)* u, . This neighborhood Q, is called the 

domain of influence of x, or circle of influence in two dimensions. The domain of 
influence of a point is shown in Fig.3.1. 



Fig. 3.1 Cell structure used for the integration for EFG and domain of 
influence around the quadrature points. 


Finding the minimum of J(x) with respect to a(x) leads to a set of linear equations 


A(x)a(x) ~ C(x)u 


(3.4) 
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For example, for linear approximation 


A (x) = J w(x- X, )p{x, )p T (x , ) 

7=1 


(3.5) 
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(3-6) 


c(x) = [w(x - X, )p(x , ), w(x - x 2 )p(x 2 ),.... w(x - x n )p{x n )] 


(3.7) 


w(x - X, ) 

■ 1 ~ 

X, 

■ w(x - x 2 ) 

■ 1 ■ 

x 2 

s: 

* 

1 

' 1 " 

x „ 



y,_ 


y 2 _ 


_y„_ 



u = \u y ,u 2 ,....u n ] 


(3.8) 


The matrix A(x) is often called the moment matrix and C(x) is a matrix of 
column vectors. So eq. 3.4 can be solved for a(x) to yield 


a(x) = A~ l (x)C(x)u 


(3.9) 


which can be substituted into eq. 3.1 to yield an approximation in terms of the 
nodal coefficients 


u ''( x ) = 1[j P T X*)^ _l ( x Kv ( x ) u i (3.10) 

7*1 

where C,(x ) is the I" 1 column of C(x)and u, is the nodal coefficient for the i"’ 
neighbor of x. Defining the shape function as f,(x) = p‘ (x^'^x^^x) 
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Allows the approximation to be written as 


u/ ’( x ) = X^/(*)"/ (3.11) 

/= i 

which is a form similar to those in the finite element method. 

It should be noted that these shape functions do not satisfy the condition 

(/>,{xj)*6 u (3.12) 

where 8„ is the Kronecker delta. In the finite element method, where is equation of 

essential boundary conditions can easily be applied, imposing essential boundary 
conditions is difficult in the EFG method. In the EFG method the displacement at a 
point Xj on the boundary node is given by 

»(*,)=!>,(*>, < 313 > 

/« 1 

where 7 = 1 to n are all nodes in the domain of influence of x, ; because of 
u(x,)*u lt so setting u, =0 on the boundary does not ensure that «(x) = 0 on the 
boundary. 

To determine the derivative form the displacement; it is necessary to obtain 
the shape function. The spatial derivatives of the shape functions, computed by the 
chain rule, are 

where A~',i =-A~ , A / A~' note that the second term in equation is expensive to 
compute because of the term A-, 1 . 


f 
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3.3 Shape Function And its Derivatives 


To compute the shape functions from the previous equation it is necessary to invert 
the A matrix. In one dimension, this operation is not computationally expensive, 
but in two or three dimensions it becomes quite difficult. An alternative approach 
that is much more computationally efficient was proposed by Belystchko al (1996). 


This method involves the LU decomposition of the A matrix. 

The shape function is equation can be written as 

0,(x)= P T ( x ) A ~ ] ( x )C,{x) 

= y T {x)C l {x) (3-15) 

with corresponding derivatives 

0,j(x) = r T j(x)C I (x)+y r {x)C, i (x) ' ( 3 - 16 ) 

Comparing the terms of eq. 3.15 leads to the relationship 

A(x)y(x) = p(x) (3-17) 


The coefficients y(x) can be obtained by an LU decomposition of A{x) and 
back substitution. This requires fewer computations than a full inversion of A(x), 
which is required to form the shape functions with equation of </>,{x). 

The derivatives of y(x) is required to compute the shape function derivatives 
in equation ^ /; (x). Taking the derivatives of eq. 3.17. 

A,, (x)y (x) + A(x)y, (x) = pXx) ( 3 • 1 

and rearranging known terms of the above equation will be, 

A {x)y l (x) = p, (*) - Aj {x)y{x) ( 3 - ' 1 9 ) 
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Using the LU decomposition of a(x), which is known from the solving 

equation, y, (x) can be computed with only a back substitution. Higher order 

derivatives can be easily obtained by repeating the procedure used for computing 
the first derivative. 

3.4 Variational Function And Discretization 

In two-dimensional, the equilibrium equations and boundary conditions are 

V-a+b = 0 in Cl (3.20) 

a ' n ~ l on r, (traction boundary conditions) 

u = u on r u (essential boundary conditions) 

In which the superposed bar denotes prescribed boundary values, and n is the unit 
normal to the domain Q . Here cr = Ds is the stress vector, D is the material 
property matrix, s = V s u is the strain vector, V v is the symmetric part of the 

gradient operator, u is the displacement vector, b is the body force vector, i and 
u are the vectors of prescribed surface tractions and displacements, respectively, 
11 is a unit normal to domain, Q, r, and r„ are the portions of boundary, r where 
tractions and displacements are prescribed, respectively. 

The variational (or weak) form for equilibrium equation can be written 

sw = V j( Su : add - ^ Su ■ bdQ. - f Su ■ tdT - SW U (u ) = ° (3. .2^ 1 ) 

The term 8W u [u ) is an additional term to account for essential boundary conditions 
in a meshless method and will be discussed in the subsequent part. 

The term 5W u {u) is necessary because, in contrast to finite-element 
methods, <j), (x, ) * 5 U , so it is not sufficient to set the nodal displacements to 
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enforce essential boundary conditions. Several forms of 5W V («) are possible. In 


this work Lagrange multipliers are used, so 

8W u = 8X • (w - u]dT + J&/ • MT (3.22) 

For linear elasticity 

s = V x u (3.23) 


which can be used to write the weak form from variational form equation in terms of 
the displacement u . 

The discrete form of variational equation for a meshless method can be 
obtained by using the approximation form equation as an approximation for u and 
8u to get the system of equations 

Ku = f m (3.24) 

The stiffness matrix K and the external force vector /“' are composed of 

the sub matrices 

K u = \B,DBjdQ (3.25) 

n 

//■ v ' = \(j), tdT+\</>,bdQ (3.26) 

r, n 

where D is the constitutive matrix given by 

l-o v 0 

£)=„ , v l-o 0 plane strain 

(1 + 41-2^) 1-2 u 

2 . 
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plane stress 


D = - 


l-v 2 


1 v 

V 1 
0 0 


0 

0 

l-v 


and B ; is a matrix of shape function derivatives 




K* 0 
o K 

<t>I,y<Pl,x 


3.5 Weight Function 

Weight function are defined to be smooth and monotonically decreasing at each 
node such that the whole domain is covered. Common shapes of weight functions 
in two dimensions are circles and rectangles (see Fig.3.2). The counterparts of 
these in three dimensions are spheres and bricks. 
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Fig. 3.2 Domain of influence in two dimensions 


An important ingredient of EFGM or the meshless methods is the weight 
function, w(x) . The choice of weight function can affect the MLS approximation of 
the u h {x) . In this analysis, the Gaussian weight function has the form, 


Gaussian: w(d,)- 


exp 


/ 

'd,) 

1 

-exp 

f 

r d : 

ml 

)1 

V 

V C J 

J 


\ 

V c . 

J J 


1-exp 


fi \ 




^ c ; 


d, * d ml 
d, > d m[ 


where d, =||x-x / || is the distance from a sampling point x to a node x, , and d mI is 

the domain of influence of a node, i.e., the area for which the weight function is 
nonzero. The variable c in the Gaussian weight is used to control the dilation of the 
weight function. It is useful to define a characteristic nodal spacing, c, , which is a 
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distance such that a node will yield a minimum set of neighbors sufficient for 
regularity of the equations used to determine the approximant. The weight function 
parameters are defined in terms of c, 


d m/ — d mnx Cj 

(3.27) 

c = ac I 

(3.28) 


where d mm and a are constants. 

For the Gaussian weight function in the equation, the parameter a is a 
dilation parameter, which controls the shape of the weight function. If a is kept 
constant while d mm is increased, the shape of the weight function will not change 
and the effective domain of influence will be smaller than the actual domain of 

d 

influence. It is necessary to have the ratio — 1 — > 4.0 to avoid poorly formed shape 

a 

functions [Belystchko,1998], In addition, a >0.5 is needed for smooth shape 
functions and derivatives. 

3.6 Evaluation of Integrals 

Computing the stiffness matrix and force vector from equation requires integration 
over the domain Q , which, in two dimensions corresponds to area integration. 
Integrating the stiffness matrix and force vector requires a numerical integration 
scheme such as Gauss quadrature, which in turn, requires a subdivision of the 
domain. Since meshless methods have no inherent subdivision of the domain like 
finite elements, it is necessary to introduce subdivisions as shown in Fig. 3. 3. The 
first one shown, which matches the domain; this technique is often called an 
element quadrature .The vertices of this background mesh are often used as the 
initial array such as the nodes at the crack tip in the model shown. 

The second integration technique, which is often called cell quadrature, uses 
a background grid of cells, which is independent of the domain. At each integration 
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point it is necessary to determine if it lies inside the domain before it is used for 
integrating equation. This technique is not widely used because it does not yield 
good accuracy along curved and angled boundaries. 

To perform the quadrature, the domain of influence of each quadrature point 
must be determined. The domain of influence of a point never extends across any 
boundaries. Thus, for a quadrature point near a crack, the domain of influence of 
x a is limited to those points x, that can be connected to x Q without interesting a 

boundary of the domain. The contributions of a quadrature point to the linear 
equations depend on the nodes in the domain of influence of the point x () . 

Domain of influence 

Properly choosing the domain of influence or nodal support is an important aspect 
of meshless methods. The size of the support should be sufficiently large so that 
the moment matrix is regular and well conditioned and that the spatial distribution 
of neighbors is fairly even. On the other hand, choosing domains of influence that 
are too large leads to a great deal of computational expense in forming the 
approximations as well as assembling the stiffness matrix. Support sizes, which are 
too large, also detract from the local character of the approximation; for problems 
involving sharp gradients, some loss of accuracy is typically noted as the effect of 
the gradient is smeared. 




3.7 Smoothing of EFG approximations 

3.7.1 Introduction 

The smoothness, which is inherent in meshless methods, has two important 
features. It provides approximation of functions and their derivatives, which are 
smooth and have the same continuity as the weight function. However, in cases 
where a discontinuity is present in either the geometry or the material, this higher 
order smoothness leads to difficulties that must be addressed in order to obtain 
good accuracy. Discontinuity occurs along an interface between two materials with 
different properties, leading to discontinuous strains across the interface. This 
situation is easily modeled using finite element where the interpellants are C 0 . 

Cordes (1996) have solved problems with multiple material is using EFG by 
treating the individual material as separate bodies and joining them together with 
Lagrange multipliers. Belytschko (1997) have developed techniques in which a so- 
called jump term is included in the trial function that is capable of representing the 
discontinuity. The magnitude of the jump becomes an unknown in the Galerkin 
discretization. 

The second class of problems, which contain discontinuities, are due to 
geometrical effects. This includes cases where a boundary of the geometry is 
nonconvex, such as a body with a hole or a crack. Because of the aforementioned 
smoothness of meshless methods, special procedures are required near 
nonconvex boundaries. There are technique to handle the nonconvex problem are 
visibility criterion, diffraction, transparency and see-through. We will use the 
diffraction method. 
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(a) Support around a crack tip 



I INTERIOR 
\ HOLE IN 
\DOMAIN 




(b) Support around a hole 


Fig. 3.4 Diffraction method for construction smooth weight function 
3.7.2 Diffraction Method 


Continuous and smooth approximants near nonconvex boundaries can be 
constructed quite easily by the diffraction method Organ (1996). The nodal support 
is wrapped around nonconvex boundaries similar to the way light diffracts around 
sharp corners. This method, which has also been called the wrap around method, 
is quite general and can be used for cracks or smooth boundaries such as interior 
holes. For example, in Fig. 3.4, where a line between the node, X,, and a sampling 
point, X , intersects a crack and the tip is within the domain of influence of the 
node. The weight function distance, d , , is modified by 



where 5 , = \\X, -X c \\, s 2 (x)=||Jf- X c \\, j 0 (x) = ||^-X / || , and X, is the node, X 
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is the sampling point, and X c is the crack tip. The parameter X is used to adjust 
the distance of the support on the opposite side of the crack. 

The spatial derivatives of the weight function are computed using the chain rule: 

dw _ dw dd, 
dx, dd, dx, 


Since ^/ dd is unchanged, all that is necessary are expressions for dd ' y 


dx 


dd, 

dX: 


= X 


^1 + S 2 ^ 


r^i + s 2 

^ j 

dx, 

l *0 J 


9s o 
dx, 


where 


ds 0 _ x, -x„ 
dx, s 0 


and 


ds 2 = X, -x a 
dx, s 2 


The diffraction method works well for general nonconvex boundaries as well. 
Consider the case shown in the line between a node and a sampling point 
intersects the boundary of a hole. The tangent point between the node and the 
nonconvex boundary is used as wrap-around point, X c and earlier equation is 

used to compute the weight function distance, d, . Note that for this implementation 
of diffraction method, the segment d 2 (x) intersects the boundary. 

One drawback of the smooth approximations is that the computed crack 
opening profile is not parabolic at the tip. The shape functions wrap around the 
crack tip and the crack in effectively shortened if the smoothing effect is too large 
or the mesh is too coarse, so the maximum stress is then not at the crack tip, but is 
shifted a small distance that depends on the mesh refinement. This effect can be 
reduced by increasing X in the diffraction method. 


34 



4.8 Closure 


The moving least square methodology is reviewed and used to construct discrete 
EFG approximation. The elastic static boundary value problem is presented along 
with its associated weak form. Approximation of the solution with EFG is presented 
and topics of nodal domains of influence and integration of weak form are 
discussed. It should be noted that the definition of a meshless method is one in 
which no predefined element connectivity exists for determining the approximant. 
Some confusion invarbly arises when the meshless approximant is used in 
Galerkin method. Integration of weak form is performed by Gauss quardrature, 
which requires some sort of integration cells. Although this detracts from the 
meshless nature of the method the background cell structure by no means 
destroys it. Smoothing of EFG approximation near nonconvex boundaries is 
presented. Without smoothing, EFG approximation near nonconvex boundaries 
such as crack tips will contain discontinuities, which extend into the domain. These 
discontinuities arise due to a sampling point grazes the boundary. The diffraction 
method smoothes EFG approximation by wrapping the nodal support a short 
distance around the point at which the discontinuity would begin. 
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CHAPTER 4 


RESULTS & DISCUSSION 

The present EFG method was applied to perform fracture mechanics 
analysis of both stationary and propagating crack. Both single (mode I) and mixed 
mode (mode I & II ) conditions were considered and five examples are presented 
here. The results are organized as follows: 

1. Convergence study 

> Near tip stress 

> Centrally located crack 

■ Small 

■ Long 

2. Crack growth analysis 

> Two symmetrical cracks 

> Two unsymmetrical cracks 

4.1 Convergence Study 

A rectangular plate with an edge crack under pure tension is considered (Fig. 4.1). 
The far field applied traction is t x 1 1 unit. The mode I stress intensity factor K x 

was evaluated by using the J integral. Due to symmetry, only half of the plate was 
analyzed. The plate is divided into 10x10 rectangular cells with their nodes 
coinciding with the meshless nodes (discussed in section 3.6) to facilitate the 
numerical integration. For all the cells, a 4x4 Gauss Quadrature was used except 
for the two around the crack tip, where 9x9 Gauss Quadrature was used to handle 
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the crack tip singularity. A plane stress condition was assumed with E = 207000 
and v = 0.3. 

The value of stress intensity factor for mode I as a function of the crack 
length is given in the Table 4.1. As can be seen, the computed stress intensity 
factor agrees with the exact values of K l for all values of crack length. The 
percentage error in stress intensity factor K l is found to be equal to 1 .229%. 

AparCfrom the crack length, the accuracy of stress intensity factor also 
depends upon the additional number of nodes around the crack tip. The effect of 
additional number of nodes and basis on stress intensity factor is compare in the 
Table 4.2. It is observed, that as the additional number of nodes increases the error 
goes down (Table 4.2) and this also shows the effect of linear and quadratic basis 
assumption. In linear basis assumption, zero additional node percentage error is 
equal to 16.01% and as go on 15 additional nodes this is reduced to 1.23%. In 
quadratic basis assumption, zero additional node percentage error is equal to 
13.02% and as go on 15 additional nodes this is reduced to 1.15%. 

Near tip stress for mode I 

The near tip stress field of a crack of size, 2a=1 unit. Containing an edge 
crack of a=0.5 unit is investigated. The meshless discretization (nodes) is shown in 
Fig. 4.1 (a.). The objective is to predict the near tip stress field by meshless method 
and to compare with corresponding LFEM field. A plane stress condition was 
assumed with the same value of E and o. For evaluation of stress, linear basis 
with 15 additional nodes is used. 

Fig. 4.2 shows the plot of radial (<r, r ) and circumferential (cr,*) stress as a 

function of rta for the range of 0.005 to 0.20, with 0 + 0°. The error between 
predicted stress values from the meshless method and the exact values is within 
2%. Fig. 4.2 shows the plot of radial (crj, circumferential (a rd ) and shear (r rg ) 
stress as a function of 0 when r/a = 0.08. The error between predicted stress 
values from the meshless method and the exact values is within 2%. 
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Centrally located straight edge crack 

> Small crack 

A centrally located small cracked plate with the following initial parameters 
and location of crack is given in Fig. 4.4 

Domain =5x3.5 cm 2 

Prack length = 0.5 cm 

> Long crack 

A centrally located long cracked plate with the following initial parameters 
location of crack is given in Fig. 4.5. 

Domain = 5x3.5 cm 2 

Crack length = 1.0 cm 

To verify the robustness of method, results were compared for different 
values of crack length. 

For both the problems, the theoretical values of K x and K n for the finite 
plate after applying the finite plate correction factor agree within the range of 0.5% 
with the values obtained from meshless method Table 4.3. 

4.2 Crack Growth Analysis 

One of the biggest benefits of the meshless methods such, as EFG is their inherent 
ability to model arbitrary crack propagation due to the absence of predefined 
element connectivity. Growing a crack consists of simply adding another line 
segment at the existing crack tip as shown in Fig. 4.3. The direction and the growth 
rate are determined by the principles of the fracture mechanics, which was 
discussed in earlier chapter. 
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First of all, a node has to be located at the crack tip. This node is actually 
placed a short distance ahead of the tip of the crack. As the crack propagates, it 
will pass directly through the current crack tip node and may pass through other 
nodes as well. In this case, two nodes replace the node, one above the crack and 
another below the crack. This is preferred to deleting the node because it 
increases the spatial resolution along the crack surface. The meshes were 
constructed by simply translating the dense pattern around the crack tip through 
the fixed mesh representing the geometry, dividing any node passed by the crack 
into two and displacing them normal to the crack surface by a distance 10' 4 w. 

Two symmetrical located 

The symmetrically located straight cracks with the following initial 
parameters and the location of crack given in the Fig. 4.6. 

Domain =5x3.5 cm 2 

After each iteration, the crack profile is given the proper crack extension in 
the right direction and the whole analysis is repeated. The stress intensity factor is 
given in Table 4.4. The corresponding crack geometries at the beginning of the 
each step are shown in the Figure (4.6-4. 9). 

1 . Following the theoretical analysis for same geometry and location of crack, 
both the crack should be propagating in a straight line. This is also what was 
observed the present analysis. 

2. As crack grows, the interaction between two cracks becomes more severe 
leading to an increase in stress intensity factor. 


Two unsymmetrical located 

The unsymmetrical located straight cracks with the following initial 
parameters and location of crack in given in figure 4.10. 
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Domain =5x3.5 cm 2 


After each iteration, the crack profile is given the proper crack extension in 
the right direction and the whole analysis is repeated. The stress intensity factor is 
given in the Table (4. 5-4. 6). The corresponding crack geometries at the beginning 
of the each step are shown in Figure (4.11-4.19). 

After each iteration, the analysis is repeated by another basis, which is 
quadratic, and, compare the results in the Figure (4.20-4.23) and Table (4. 5-4.6). 

1. As the cracks grow, the interaction between the two cracks increases 
leading to an increase in the observed stress intensity factor values. 

2. The cracks tend to turn towards each other. So the intuitive reasoning that 
near the boundary, the crack will tend to grow towards the region of high 
stress concentration. 

3. The mixed mode stress intensity factor is mainly K, dominated. Mode II is 
important for the point of view of orientation of the crack growth direction. 
The crack turns primarily due to K„ . 

4. Difference between the different basis is evaluation is very less and so 
quadratic or higher degree basis can be sweep by raising the number of 
node around the crack tip. 
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Table 4.1: Stress intensity factor for the edge crack problem 


Crack length (a) 


K 

^ exact 

% Error 

0.20 

1.078 

1.086 

0.736 

0.24 

1.249 

1.278 

2.269 

0.28 

1.481 

1.493 

0.803 

0.32 

/ 

1.722 

1.737 

0.863 

0.36 

2.073 

2.021 

2.572 

0.40 

2.329 

2.358 

1.229 


Table 4.2: Stress intensity factor comparison for edge crack 


Additional 

Basis 

Node 

Linear 

Quadratic 


K. 

^*1 exact 

K, 

Kj / K exac( 

0 

1.9807 

0.8399 

2.051 

0.8698 

5 

2.2122 

0.9381 

2.284 

0.9686 

15 

2.329 

0.9877 

2.3308 

0.9885 


Table 4.3: Comparison of the exact and calculated value of S.I.F. 


Type of crack 

K, 

K 

exact 

K„ 

K' 

^ exact 

Small 

0.8955 

0.8963 ' 

-0.00072 

0.0000 

Long 

1.2757 

1.2862 

-0.00264 

0.0000 
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Table 4.4: Comparison the S.I.F. for symmetrical crack problem 


Step 

Linear 

Quadratic 


K, 

K„ 

e 

K« 


6 

1 

0.8487 

-0.00023 

0.031 

0.8503 

-0.00024 

0.032 

0.8653 

-0.00024 

0.032 

0.8679 

-0.00026 

0.034 

2 

0.8429 

-0.00017 

0.026 

0.8511 

-0.00019 

0.026 

0.8726 

-0.00048 

0.062 

0.8806 

-0.00041 

0.054 

3 

0.9216 

-0.00053 

0.068 

0.9249 

-0.00067 

0.083 

1.005 

-0.00139 

0.158 

1.0092 

-0.00147 

0.167 


Table 4.5: Comparison the S.I.F. for unsymmetrical crack problem 


Step 

Linear 

Quadratic 

% Difference 

K, 

K„ 

K, 

K„ 

K, 

Ku 

1 

0.8343 

0.0121 

0.8349 

0.0122 

-0.07192 

-0.82645 

0.8527 

0.0113 

0.8561 

0.0119 

-0.39873 

-5.30973 

2 

0.8318 

0.00561 

0.8321 

0.00567 

-0.03607 

-1.06952 

0.8588 

0.00731 

0.8602 

0.00789 

-0.16302 

-7.93434 

3 

1.0021 

-0.00162 

1.019 

-0.00145 

-1.68646 

10.49383 

1.0879 

0.0072 

1.0901 

0.0064 

-0.20222 

11.11111 

4 

1.3614 

-0.0219 

1.3671 

-0.0211 

-0.41869 

3.652968 

1 .4608 

-0.0318 

1.4691 

-0.0284 

-0.56818 

10.69182 

5 

1.7288 

-0.0397 

1.7341 

-0.0403 

-0.30657 

-1.51134 

1.8547 

-0.198 

1.8613 

-0.2028 

-0.35585 

-2.42424 

6 

2.0304 

-0.0527 

2.0706 

-0.0556 

-1.97991 

-5.50285 

2.1568 

-0.3365 

2.1812 

-0.3398 

-1.13131 

-0.98068 

7 

2.4817 

-0.2238 

2.4893 

-0.2312 

-0.30624 

-3.30652 

2.5665 

-.3191 

2.5688 

-0.3092 

-0.08962 

3.102476 

8 

2.9423 

-0.3133 

2.9512 

-0.3183 

-0.30248 

-1.59591 

2.7168 

-0.3087 

2.7218 

-0.3155 

-0.18404 

-2.20279 

9 

3.2987 

-0.4210 

3.3012 

-0.4297 

-0.07579 

-2.06651 

2.992 

-0.3942 

3.002 

-0.4011 

-0.33422 

-1.75038 

10 

3.6791 

-0.4987 

3.6912 

-0.5135 

-0.32888 

-2.96772 

3.2147 

-0.3782 

3.3128 

-0.3841 

-3.05161 

-1.56002 
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Table 4.6: Comparison the direction for unsymmetrical crack 
problem 


Step 

Linear 

Quadratic 

% Difference 

e 

9 

e 

1 

-1.66 

- 1.674 

-0.84337 

-1.52 

-1.592 

-4.73684 

2 

-0.777 

-0.781 

-0.5148 

-0.975 

-1.051 

-7.79487 

3 

0.185 

0.163 

11.89189 

' -0.756 

-0.672 

11.11111 

4 

1.845 

1.768 

4.173442 

2.49 

2.213 

11.1245 

5 

2.52 

2.66 

-5.55556 

11.92 

12.156 

-1.97987 

6 

2.99 

3.0702 

-2.68227 

16.96 

16.936 

0.141509 

WBM 

10.15 

10.437 

-2.82759 

13.77 

13.355 

3.013798 

8 

11.89 

12.04 

-1.26156 ' 

12.64 

12.89 

-1.97785 

9 

14.10 

14.366 

-1.88652 

14.52 

14.719 

-1.37052 

■BMW 

14.92 

15.275 

-2.37936 

13.07 

12.893 

1.354246 
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Fig.4.1 Single edge cracked plate 
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Fig. 4.2(a.) Near crack tip stress field for <9=0 
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Fig. 4.2(d.) angular variation of stress for r/a = 0.08 



mesh arrangement 


Fig. 4.3 Meshes for the progressive growth of edge crack problem 
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Fig. 4.4 Straight edge small crack 



Fig. 4.5 Straight edge long crack 
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Step 1 



Fig 4.6 Problem 1 growth step 1 


Step 2 


Fig 4.7 Problem 1 growth step 2 
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Step 3 



Step 4 



Fig 4.8 Problem 1 growth step 3 



Fig 4.9 Problem 1 growth step 4 
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Stepl 




Step2 



Fig 4.10 Problem 2 growth step 1 



Fig 4.11 Problem 2 growth step 2 
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Step 3 



Fig 4.12 Problem 2 growth step 3 


Step 4 



Fig 4.13 Problem 2 growth step 4 
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Step 5 



Step 6 


Fig 4.14 Problem 2 growth step 5 
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Step 7 



Step 8 


Fig 4.16 Problem 2 growth step 7 
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Fig 4.17 Problem 2 growth step 8 
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Step 10 


Fig 4.18 Problem 2 growth step 9 


Fig 4.19 Problem 2 growth step 10 



comparsion 



Fig. 4.21 Problem 2 comparison of angle at different basis 
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Fig. 4.22 Problem 2 comparison of S.I.F. at different basis 
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Fig. 4.23 Problem 2 comparison of angle at different basis 
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CHAPTER 5 


CONCLUSIONS 


5.1 Conclusions 

An Element Free Galerkin method (EFG) has been developed to analyze 2-D 
fracture problems. This method can determine crack singularity in mode I & II and 
predict crack propagation. Several studies have been presented of different crack 
geometries and EFG method parameter. 

Based on the studies the following conclusions are drawn: 

1. The application EFG to fracture is very accurate and computationally very 
efficient. The SIF (K,) is accurate and has error below 2% for linear basis and 1% 
for quadratic basis assumption. 

2. The crack propagation direction is accurate and has error less than 2%. This 
became possible with moderate mesh configuration. 

3. Unlike BEM or FEM, the EFG method does not need the crack to propagate only 
along the element edges. That way, any arbitrary propagation direction of the crack 
can be predicted. 

5.2 Suggestions For The Future Work 

1 . Further studies, regarding the other basis functions and weighting functions can 
be explored for more efficient algorithms. 

2. The method can be extended to analyze dynamic fracture problems and elasto- 
plastic material fracture problems. 
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